// -*- C++ -*------------------------------------------------------------------
// $Id: PdbLoader.cc,v 1.22 2008-08-03 21:26:35 biocomp Exp $
//
//  Class:             PdbLoader
//
//  Base Class(es):    Loader
//
//  Derived Class(es): -
//
//  Containing:        -
//
//  Author:            Silvio Tosatto
//
//  Project Name:      Victor
//
//  Date:              09/99
//
//  Reviewed By:       <Name>
//
//  Description:
//    Loads components (Atoms, Groups, etc.) in standard PDB format.
//    NB: Only loadSpacer() is currently implemented.
//
// ---------------------------------------------------------------------------

// Includes:
#include <PdbLoader.h>
#include <IoTools.h>
#include <vector3.h>
#include <AtomCode.h>
#include <AminoAcid.h>
#include <String2Number.h>

// Global constants, typedefs, etc. (to avoid):

using namespace Biopool;

// CONSTRUCTORS/DESTRUCTOR:

// PREDICATES:

// -*- C++ -*-----------------------------------------------------------------
//
//  Method:        PdbLoader::getMaxModels
//
//  Author:        Silvio Tosatto 
//
//  Date:          08/05
//
//  Description:
//    Reads in the maximum allowed number of NMR models, zero otherwise.
//
// ----------------------------------------------------------------------------
unsigned int 
PdbLoader::getMaxModels()
{
  input.clear();  // reset file to previous content 
  input.seekg(0);

  unsigned int max = 0;

  while (input)
    {
      string atomLine = readLine(input);
      string tag = atomLine.substr(0,6);

      if (tag == "MODEL ")
	max++;
    }

  return max;
}


// -*- C++ -*-----------------------------------------------------------------
//
//  Method:        PdbLoader::getMaxModelsFast
//
//  Author:        Luca Stoppa 
//
//  Date:          08/05
//
//  Description:
//    Reads in the maximum allowed number of NMR models, zero otherwise.
//    This is a faster version because we only need to read the line with
//    EXPDTA
//
// ----------------------------------------------------------------------------
unsigned int 
PdbLoader::getMaxModelsFast()
{
	input.clear();
	input.seekg( 0 );
	while ( input )
	{
		string line=readLine( input );
		if ( line.substr(0,6)=="EXPDTA" )
		{
			istringstream str( line );
			string header, nmr, structures;
			unsigned int models;
			str >> header >> nmr >> models >> structures;
			return models;
		}
	}
	
	return 0;
}


// -*- C++ -*-----------------------------------------------------------------
//
//  Method:        PdbLoader::getAllChains
//
//  Author:        Silvio Tosatto 
//
//  Date:          08/05
//
//  Description:
//    Returns all available chain IDs for a PDB file.
//
// ----------------------------------------------------------------------------
vector<char> 
PdbLoader::getAllChains()
{
  vector<char> res;
  char lastChain = ' ';
  input.clear();  // reset file to previous content 
  input.seekg(0);
  unsigned int max = 0;

  while (input)
    {
      string atomLine = readLine(input);

      if (atomLine.substr(0,5) == "MODEL")
	max++;
      if (max > 1) // only consider first model: others duplicate chainIDs
	break;

      // check for new chains containing amino acids
      if ((atomLine.substr(0,4) == "ATOM") 
	  && (atomLine.substr(21,1)[0] != lastChain)
	  && (aminoAcidThreeLetterTranslator(atomLine.substr(17,3)) != XXX))
	{
	  lastChain = atomLine.substr(21,1)[0];
	  res.push_back(lastChain);
	}
    }

  return res;
}



// MODIFIERS:

// -*- C++ -*-----------------------------------------------------------------
//
//  Method:        PdbLoader::loadLigandSet
//
//  Author:        Silvio Tosatto 
//
//  Date:          02/02
//
//  Description:
//    Loads a LigandSet in PDB format.
//    Attention: May not load all structures correctly.
//
// ----------------------------------------------------------------------------
void 
PdbLoader::loadLigandSet(LigandSet& ls)
{
  input.clear();  // reset file to previous content 
  input.seekg(0);

  PRINT_NAME;

  while (!checkForKeyword(input, "HETATM"))
      if (!input)
	  return;
      else
	  skipToNewLine(input);

  int oldAaNum = -2000;
  Ligand* lig = NULL;
  bool continueReading = true;

  do {  // read all entries

    int tNum, aaNum;
    string type, aaCode, tmpType;
    vgVector3<double> coord;
    input >> tNum >> type;
    if (type.length() > 4)
      {
	aaCode = type.substr(3);
	type.erase(3);
      }
    else
      input >> aaCode;          // read a single entry
    if (aaCode.length() > 3)
      {
	aaCode = aaCode.substr(1);
      }

    if (noWater && (aaCode == "HOH"))   // skip water if requested
    {
	skipToNewLine(input);
	continueReading = checkForKeyword(input, "HETATM");
	if ((!continueReading) && (input))  // skip single "ANISOU" entries
	{
	    skipToNewLine(input);  
	    continueReading = checkForKeyword(input, "HETATM");
	}
	continue;
    }

    if (readNumber(input, aaNum))
      input >> coord; 
    else
      {
	if (!permissive)
	  { ERROR("File is not in proper PDB format for ligands?!", exception); }
	else
	  {
	    valid = false;
	    return;
	  }
      }

    // read B-factor:
    double /*occ,*/ Bfac = 0.0; 
    string rest;
//    input >> occ >> Bfac; 
    rest = readLine(input);

    if (rest[4] != ' ')
	Bfac = stod(rest.substr(3,7));
    else
	Bfac = stod(rest.substr(4,6));
	
//      skipToNewLine(input);  

    if (aaNum != oldAaNum) // check if a new ligand has to be created
      {
	lig = new Ligand;
	lig->setType(aaCode);

	ls.insertComponent(lig);
//  	if (oldAaNum == -2000)
//  	  ls.setStartOffset(aaNum);
	oldAaNum = aaNum;
      }
    Atom at;
    at.setType(type);
    at.setCoords(coord);
    at.setBFac(Bfac);
//      if ( (!noHAtoms) || (isHeavyAtom(at.getCode())) )
	lig->addAtom(at);

    eatComment(input);
    continueReading = checkForKeyword(input, "HETATM");
    if ((!continueReading) && (input))  // skip single "ANISOU" entries ???
    {                             // not sure if this case exists for ligands  
	skipToNewLine(input);  
	continueReading = checkForKeyword(input, "HETATM");
    }
  } while (continueReading);

}



// -*- C++ -*-----------------------------------------------------------------
//
//  Method:        PdbLoader::loadSpacer
//
//  Author:        Silvio Tosatto 
//
//  Date:          01/03
//
//  Description:
//    Loads a spacer in PDB format.  (ver 2.0)
//    Completely re-written to be as fault-tolerant as possible.
//    It is now possible to load structures with multiple chains, residues 
//    and atoms. Alternate configurations are largely ignored.
//  Note: We changed the order of functions calls inverting seekg() & clear()
// ----------------------------------------------------------------------------
void 
PdbLoader::loadSpacer(Spacer& sp)
{
	// Pay attention to the fact that
	// if you execute a 
	input.clear();  // reset file to previous content 
	input.seekg( 0, ios::beg );

	PRINT_NAME;
	if (checkForKeyword(input, "HEADER"))
    {
		//cout <<"yyyyy" <<endl;
		string name;   // title line
		name = readLine(input);
		sp.setType(name);
    }
	// -----
	
	string atomLine;
	bool firstAtomLineRead = false;

	if (!noSecondary)
    {
		helixCode = "";
		sheetCode = "";
		unsigned int start, end;
		
		while (!input.eof())
		{
			atomLine = readLine(input);
			string tag = atomLine.substr(0,6);

			if ( (tag == "MODEL ") || (tag == "ATOM  "))
			{
				firstAtomLineRead = true;
				//cout <<"I did it" <<endl;
				break;
			}

			if (tag == "HELIX ")
			{ // read helix entry:
	      //****************************************
//  	      cout << "AtomLine:  ::" << tag << "::"
//  		   << atomLine.substr(19,1) << "::" 
//  		   << atomLine.substr(21,4) << "::" 
//  		   << atomLine.substr(33,4) << "::" << endl;
	      //****************************************
				start = stoui(atomLine.substr(21,4));
				end = stoui(atomLine.substr(33,4));
	      //****************************************
//  	      cout << "AtomLine:  ::" << tag << "::"
//  		   << start << "::" << end << "::" << endl;
	      //****************************************
				helixData.push_back(pair<const unsigned int, unsigned int>(start, end));
				helixCode += atomLine.substr(19,1).c_str()[0];
			}
			else if (tag == "SHEET ") 
			{ // read sheet entry:
	      //****************************************
//  	      cout << "AtomLine:  ::" << tag << "::"
//  		   << atomLine.substr(21,1) << "::" 
//  		   << atomLine.substr(22,4) << "::" 
//  		   << atomLine.substr(33,4) << "::" << endl;
	      //****************************************
				start = stoui(atomLine.substr(22,4));
				end = stoui(atomLine.substr(33,4));
	      //****************************************
//  	      cout << "AtomLine:  ::" << tag << "::"
//  		   << start << "::" << end << "::" << endl;
	      //****************************************
				sheetData.push_back(pair<const unsigned int, unsigned int>(
				   start, end));
				sheetCode += atomLine.substr(21,1).c_str()[0];
			}
		} // end-while
    }// end-nosecondary

  int oldAaNum = -2000;
  int oldWarnAa = -2000;
  int aaNum = -2000;

  char chainLoad = chain;

  AminoAcid* aa = NULL;
  bool continueReading = true;
  bool notModelWarned = true;

  bool inAtoms = false;
	//input.seekg( 0, ios::beg );
	//input.clear();  // reset file to previous content 
  unsigned int readingModel = model;


  do {  // read all entries
	
    if (firstAtomLineRead)
      firstAtomLineRead = false;
    else
      atomLine = readLine(input);
    //****************************************
//      cout << "AtomLine:  ::" << atomLine << "::" << endl;
//    	 << atomLine.substr(6,5) << "::" << endl;
    //****************************************

    string tag = atomLine.substr(0,6);
    if (tag == "MODEL ")
      {
	if (notModelWarned)
	  {
	    if (model == 999)
	      model = stoui(atomLine.substr(6,10));
	    else
	      if (verbose)
		{
		  cout << "Loading model " << model << "\n";
		}

	    notModelWarned = false;
	  }

	readingModel = stoui(atomLine.substr(6,10));
	continue;
      }

//      if (tag == "ENDMDL")
//        {
//  	continueReading = false;
//  	continue;
//        }

    if ((tag != "ATOM  ") || (readingModel != model))
      continue;

    int atNum = stoi(atomLine.substr(6,5));

    // adjust atType to match just non-empty characters:
    int i1 = 0;
    for (int i = 0; i < 5; i++)
      if (atomLine[11+i] == ' ')
	i1 = i+1;
      else
	break;
    string atType = atomLine.substr(11+i1,5-i1);
    while (atType[atType.size()-1] == ' ')
      atType.erase(atType.size()-1);
    
    char altAtID = atomLine.substr(16,1)[0];
    string aaType = atomLine.substr(17,3);
    char chainID = atomLine.substr(21,1)[0];
    aaNum = stoi(atomLine.substr(22,4));
    char altAaID = atomLine.substr(26,1)[0];

    vgVector3<double> coord;
    coord.x = stod(atomLine.substr(30,8));
    coord.y = stod(atomLine.substr(38,8));
    coord.z = stod(atomLine.substr(46,8));

    double bfac = 0.0;
    if (atomLine.length() >= 66)
      bfac = stod(atomLine.substr(60,6));

    //****************************************
//      double occ = stod(atomLine.substr(55,5));
//      cout << "AtomLine:  ::" << tag << "::" 
//  	 << atNum << "::" << atType << "::" << altAtID << "::"
//  	 << aaType << "::" << chainID << "::" << aaNum << "::" 
//  	 << altAaID << "::" 
//  	 << coord.x << "::" << coord.y << "::" << coord.z << "::" 
//  	 << occ << "::" << bfac << "::" 
//  	 << endl;
    //****************************************

    inAtoms = true;

    if ((chainLoad == ' ') && (chainID != ' '))
      {
	if (verbose)
	  {
	    cout << "Loading chain " << chainID << "\n";
	  }
	chainLoad = chainID;
	chain = chainID;
      }

    if (altAaID != ' ')
      {
	if (oldWarnAa != aaNum)
	  {
	    if (verbose)
	      {
		cout << "Warning: Skipping extraneous amino acid entry " 
		     << aaNum << altAaID << ".\n"; 
	      }
	    oldWarnAa = aaNum;
	  }
	continue;     // skip alternate amino acid entries 
      }

    if ((altAtID != ' ') && (altAtID != altAtom))
      continue;     // skip alternate atom entries if not correct
                    // (usually only includes 'A' atoms)

    if (chainID != chainLoad)
      continue;

    if (aaType == "ACE") // skip N-terminal ACE groups 
      continue;

    if (atType == "D")
      {                // take care of deuterium atoms
	cerr << "--> " << atType << "\n";
	atType = "H";
      }

    if (aaNum != oldAaNum) // check if a new aminoacid has to be created
      {
	aa = new AminoAcid;
	aa->setType(aaType);
	aa->getSideChain().setType(aaType);
	if (!noConnection)
	    aa->getSideChain().setBackboneRef(aa);
	sp.insertComponent(aa);
	if (oldAaNum == -2000)
	  sp.setStartOffset(aaNum);
	else
	  {  // track missing AAs ('gaps') in PDB file
	    for (int i = oldAaNum+1; i < aaNum; i++)
	      sp.addGap(i);
	  }
	oldAaNum = aaNum;
      }
    Atom at;

    at.setNumber(atNum);

    at.setType(atType);
    at.setCoords(coord);
    at.setBFac(bfac);
    if ( (!noHAtoms) || (isHeavyAtom(at.getCode())) )
      {
	if (inSideChain(*aa, at))
	  aa->getSideChain().addAtom(at);
	else
	  aa->addAtom(at);
      }

  } while ((continueReading) && (input));

//******************************
//   cout << "ending main loop..." << endl;
//******************************


  if (sp.sizeAmino() == 0)
    {
      if (permissive)
	{
	  valid = false;
	  return;
	}
      else
	{
	  if (!inAtoms)
	    { ERROR("Invalid file format. No ATOM entry found.", exception); }
	  else
	    { ERROR("No ATOM entry with matching chain ID found.", exception); }
	}
    }

  unsigned int lastDel = 0;

  // correct ''fuzzy'' (i.e. incomplete) residues:
  for (unsigned int i = 0; i < sp.sizeAmino(); i++)
    {
      if ( (!sp.getAmino(i).isMember(O)) ||
	   (!sp.getAmino(i).isMember(C)) ||
	   (!sp.getAmino(i).isMember(CA)) ||
	   (!sp.getAmino(i).isMember(N)) )
	{

	  // remove residue
	  sp.deleteComponent(&sp.getAmino(i));

	  // increment startoffset if first residue in chain
	  if (i == lastDel)
	    {
	      sp.setStartOffset(sp.getStartOffset()+1);
	      lastDel++;
	    }

	  if (verbose)
	    {
	      cout << "Warning: Residue number " 
		   << sp.getPdbNumberFromIndex(i) 
		   << " is incomplete and had to be removed.\n";
	    }
 	}
   }


  if ((!noConnection) && (sp.sizeAmino()))
      if (!setBonds(sp))  // connect atoms...
	{
	  valid = false;
	  return;
	} 

  if (sp.sizeAmino()) // correct position of leading N atom
    {
      sp.setTrans(sp.getAmino(0)[N].getTrans());
      vgVector3<double> tmp(0.0, 0.0, 0.0);
      sp.getAmino(0)[N].setTrans(tmp);
      
      sp.getAmino(0).adjustLeadingN();
    }

  assignSecondary(sp);

}


// -*- C++ -*-----------------------------------------------------------------
//
//  Method:        PdbLoader::setBonds
//
//  Author:        Silvio Tosatto 
//
//  Date:          09/99
//
//  Description:
//    Private helper function to set bond structure after loading the spacer.
//
// ----------------------------------------------------------------------------
bool
PdbLoader::setBonds(Spacer& sp)
{
  sp.getAmino(0).setBondsFromPdbCode(true);
  for (unsigned int i = 1; i < sp.size(); i++)
    if (!sp.getAmino(i).setBondsFromPdbCode(true, &(sp.getAmino(i-1)) ))
      return false;
  return true;
}


// -*- C++ -*-----------------------------------------------------------------
//
//  Method:        PdbLoader::inSideChain
//
//  Author:        Silvio Tosatto 
//
//  Date:          09/99
//
//  Description:
//    Private helper function to determine if atom is backbone or sidechain. 
//
// ----------------------------------------------------------------------------
bool 
PdbLoader::inSideChain(const AminoAcid& aa, const Atom& at)
{
  if (isBackboneAtom(at.getCode()))  
    return false;
  if ( (at.getType() == "H") || (at.getType() == "HN") 
       || ((at.getType() == "HA") && (!aa.isMember(HA)))
       || (at.getType() == "1HA") || (at.getType() == "1H")
       || (at.getType() == "2H") || (at.getType() == "3H") )
    return false;  // special case for GLY H (code HA)
  return true;  // rest of aminoacid is its sidechain
}


// -*- C++ -*-----------------------------------------------------------------
//
//  Method:        PdbLoader::assignSecondary
//
//  Author:        Silvio Tosatto 
//
//  Date:          05/00
//
//  Description:
//    Assigns the secondary structure (if any) loaded in loadSecondary() or
//    uses AminoAcid's setStateFromTorsionAngles() if no structure was found.
//
// ----------------------------------------------------------------------------

void
PdbLoader::assignSecondary(Spacer& sp)
{
  if (helixData.size() + sheetData.size() == 0)
    {
      sp.setStateFromTorsionAngles();
      return;
    }

  for (unsigned int i = 0; i < helixData.size(); i++)
    if (helixCode[i] == chain)
      for (unsigned int j = helixData[i].first; 
	   j <= const_cast<unsigned int&>(helixData[i].second); j++)
	{ // important: keep ifs separated to avoid errors
	  if (j < sp.maxPdbNumber()) 
	    if (!sp.isGap(sp.getIndexFromPdbNumber(j)))
	    sp.getAmino(sp.getIndexFromPdbNumber(j)).setState(HELIX);
	}
  for (unsigned int i = 0; i < sheetData.size(); i++)
    if (sheetCode[i] == chain)
      for (unsigned int j = sheetData[i].first; 
	   j <= const_cast<unsigned int&>(sheetData[i].second); j++)
	{ // important: keep ifs separated to avoid errors
	if (j < sp.maxPdbNumber())
	  if (!sp.isGap(sp.getIndexFromPdbNumber(j))) 
	      sp.getAmino(sp.getIndexFromPdbNumber(j)).setState(STRAND);
	}
}
